Digital Filtering - A Fundamental Pre-Processing Step
Difficulty Level:
Tags pre-processā˜filter

The acquired electrophysiological signals have always two intrinsic components. The signal we really want to acquire/study and noise, i.e. the acquisition component that is not relevant for the study purposes.

Noise can have different origins, such as in random events or due to voluntary/involuntary movements of the subject under analysis that affect the signal acquisition .

So, filtering is a fundamental step that needs to be applied to the signal, in order to ensure the maximisation of signal to noise ratio . Filtering can be achieved by hardware, having the analogical systems a great relevance, or by software using digital filters.

In this Jupyter Notebook it will be demonstrated how to digital filter the signal.

1 - Importation of the needed packages

InĀ [1]: 
# biosignalsnotebooks own package for loading and plotting the acquired data
import biosignalsnotebooks as bsnb

# Scientific package
import numpy

# Signal Processing package where digital filtering functions are included
from novainstrumentation.tools import plotfft
from novainstrumentation.filter import bandpass
InĀ [2]: 
# Hide warnings.
# https://groups.google.com/a/continuum.io/forum/#!topic/bokeh/h817iNS2twk
from IPython.display import HTML
HTML('''<script>
code_show_err=false; 
function code_toggle_err() {
 if (code_show_err){
 $('div.output_stderr').hide();
 } else {
 $('div.output_stderr').show();
 }
 code_show_err = !code_show_err
} 
$( document ).ready(code_toggle_err);
</script>
To toggle on/off output_stderr, click <a href="javascript:code_toggle_err()">here</a>.''')
Out[2]:
To toggle on/off output_stderr, click here .

2 - Load of acquired ECG data

InĀ [3]: 
# Load of data
data, header = bsnb.load_signal("ecg_4000_Hz", get_header=True)

In the following cell, some relevant information is stored inside variables. This relevant information includes the mac-address of the device, channel number and signal acquisition parameters such as resolution and sampling rate.

For a detailed explanation of how to access this info, the "Signal Loading - Working with File Header" Notebook should be consulted.

InĀ [4]: 
mac = "00:07:80:D8:A7:F9" # Mac-address
ch = "CH1" # Channel
sr = header[mac]["sampling rate"] # Sampling rate
resolution = header[mac]["resolution"] # Resolution (number of available bits)

3 - Generation of signal power spectrum by Fast Fourier Transform (FFT)

With this step is possible to observe the frequency composition of ECG signal.

InĀ [5]: 
# Acquired data
signal = data[mac][ch]

# Power spectrum
freq_axis_1, power_spect_1 = plotfft(signal, sr)

4 - The informational content of ECG signal is typically contained inside the frequency band [0.5, 40] Hz

With the next representation we can conclude that exist some unwanted information out of this frequency band.

InĀ [6]: 
fig_1 = bsnb.plot_informational_band(freq_axis_1, power_spect_1, signal, sr, band_begin=0.5, band_end=40, legend="ECG Power Spectrum", 
                                     x_lim=[0, 100], y_lim=[0, 5e6])

5 - Application of a pass-band filter in order to be excluded the unwanted information out of [0.5, 40] Hz frequency band

InĀ [7]: 
# Digital filtering with a low cutoff frequency f1 of 0.5 Hz and high cutoff frequency f2 of 40 Hz
filter_signal_1 = bandpass(signal, f1=0.5, f2=40, order=1, fs=sr, use_filtfilt=True)

# Power spectrum
freq_axis_2, power_spect_2 = plotfft(filter_signal_1, sr)

6 - Comparison of the power spectrum of original and filtered signal

InĀ [8]: 
bsnb.plot_before_after_filter(signal, sr, band_begin=0.5, band_end=40, x_lim=[0, 100], y_lim=[0, 5e6])
Out[8]: 
[[Figure(id='27a6349c-f596-4c36-b603-7984169f9c48', ...),
  Figure(id='37ccfc7f-6f08-41ca-abd7-9aeb325b8961', ...)]]

In this first filtering attempt we used a first order filter (input argument order=1). It can be seen, in the previous figure, that some unwanted information have been removed, unfortunately no filter has an ideal behavior, so despite we specify a high cutoff frequency of 40 Hz, some information above this threshold is maintained after filtering.

The good news are that components greater than 80 Hz are almost completely removed.

The filter performance can be improved by increasing the filter order, because the higher the filter order is, more quickly the transition between the pass and stop band will be. The transition band will be smaller because of a higher attenuation rate (-20 x order dB/decade).

However, the filter order must be chosen with precaution in order to avoid system instability.

InĀ [9]: 
bsnb.plot_low_pass_filter_response()

7 - Repetition of the filtering stage but using a higher filter order

InĀ [10]: 
# Digital filtering with a low cutoff frequency f1 of 0.5 Hz and high cutoff frequency f2 of 40 Hz
filter_signal_2 = bandpass(signal, f1=0.5, f2=40, order=3, fs=sr, use_filtfilt=True)

# Power spectrum
freq_axis_3, power_spect_3 = plotfft(filter_signal_2, sr)
InĀ [11]: 
bsnb.plot_before_after_filter(signal, sr, band_begin=0.5, band_end=40, order=3, x_lim=[0, 100], y_lim=[0, 5e6], 
                              horizontal=False)
Out[11]: 
[[Figure(id='bfef6fbe-2c71-4dc0-ac58-f3030e7b4256', ...)],
 [Figure(id='49d648e1-94cd-481f-b7fa-76b360c0101b', ...)],
 [Figure(id='fe65229b-59b1-4a8c-a398-c58d68655448', ...)],
 [Figure(id='82966545-7d2e-4a5b-b4aa-4828d638d17f', ...)]]

When the noise level is low, it may be difficult to observe its influence in time domain. In order to the digital filtering stage produce a visual effect in time domain, the noise level needs to be high.

So we will add some artificial noise do the signal and see the great impact of digital filtering.

E1 - Addition of artificial noise

InĀ [12]: 
# Noise samples and translation of the baseline 
baseline = numpy.average(signal)
baseline_shift = 0.50 * baseline
noisy_signal = signal + numpy.random.normal(0, 1000, len(signal)) + baseline_shift

E2 - Noisy signal representation

InĀ [13]: 
# Plotting of power spectrum    
bsnb.plot(numpy.linspace(0, len(noisy_signal) - 1, len(noisy_signal)), noisy_signal, x_axis_label='Sample Number', 
          y_axis_label='Raw Data', title="Noisy Signal", y_range=(-1e4, 6e4))

E3 - Digital Filtering Stage

InĀ [14]: 
# Digital filtering with a low cutoff frequency f1 of 0.5 Hz and high cutoff frequency f2 of 40 Hz
noisy_signal_filter = bandpass(noisy_signal, f1=0.5, f2=40, order=3, fs=sr, use_filtfilt=True)

E4 - Comparison of noisy and filtered signal in time domain

InĀ [15]: 
bsnb.plot([numpy.linspace(0, len(noisy_signal) - 1, len(noisy_signal))]*2, [noisy_signal, noisy_signal_filter],
          grid_plot=True, grid_lines=1, grid_columns=2, x_axis_label='Sample Number', y_axis_label='Raw Data', 
          title=["Noisy Signal", "Filtered Signal"], y_range=(-1e4, 6e4))

**Auxiliary Code Segment (should not be replicated by the user)**

InĀ [16]: 
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()

.................... CSS Style Applied to Jupyter Notebook .........................
Out[16]: